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ABSTRACT 

A kilohertz quasi-periodic oscillation (kHz QPO) is an observationally robust high-frequency timing 
feature detected from neutron star low- mass X-ray binaries (LMXBs). This feature can be very useful 
to probe the superdense core matter of neutron stars, and the strong gravity regime. However, although 
many models exist in the literature, the physical origin of kHz QPO is not known, and hence this 
feature cannot be used as a tool yet. The energy dependence of kHz QPO fractional rms amplitude is 
an important piece of the jigsaw puzzle to understand the physical origin of this timing feature. It is 
known that the fractional rms amplitude increases with energy at lower energies. At higher energies, the 
amplitude is usually believed to saturate, although this is not established. We combine tens of lower kHz 
QPOs from a neutron star LMXB 4U 1728-34 in order to improve the signal-to-noise-ratio. Consequently, 
we, for the first time to the best of our knowledge, find a significant and systematic decrease of the 
fractional rms amplitude with energy at higher photon energies. Assuming an energy spectrum model, 
blackbody+powerlaw, we explore if the sinusoidal variation of a single spectral parameter can reproduce 
the above mentioned fractional rms amplitude behavior. Our analysis suggests that the oscillation of 
any single blackbody parameter is favored over the oscillation of any single powerlaw parameter, in order 
to explain the measured amplitude behavior. We also find that the quality factor of a lower kHz QPO 
does not plausibly depend on photon energy. 

Subject headings: accretion, accretion disks — methods: data analysis — stars: neutron — X-rays: 
binaries — X-rays: individual: 4U 1728-34 — X-rays: stars 



1. INTRODUCTION 

Kilohertz quasi-periodic oscillations (kHz QPOs), dis- 
covered in 1996 with Ross i X-ray Timing Explorer {R XTE; 
Ivan der Klis et al. I (|1996D ; IStrohmaver et al. I(|1996h ). are 
the fastest variability features known till date in low- 
mass X-ray binaries (LMXBs). They have been observed 
from many neutron star LMXBs, although for a given 
source they are no t always detected (jBhattacharvva I feOlO; 
van der Klisll2006h . These QPOs often occur in pairs, and 
the twin peaks usually move together in the frequency 
range of ~ 400 — 1200 Hz. The higher frequency QPO 
is known as the upper kHz QPO (frequency u u ), and 
the lower frequency QPO is called the lower kHz QPO 
(frequency v{). The frequency differences of these QPOs 
(Av = v u — v\) tend to cl uster around the ne utron star 
spin frequency or half of it (Ivan der Klisl (12006 ) an d refer- 
ences therein; but see iMendez and Bellonil (|2007l )). The 
high frequencies of kHz QPOs point toward the dynami- 
cal timesc ale within a few Schw a rzschild radii of the neu- 
tron star (iBarret et al. 1 l2005aHbl. 120061: Ivan der Klisll2006l: 
iBhattacharvva II2010D . Therefore, this observationally ro- 
bust timing feature can be useful (1) to measure the 
neutron star parameters, which, in turn, is essential to 
unde rstand the nature of s uper-dense degenerate matter 
' ;.g., Bhattacharyya I (|2010l ) and references therein; iPetri I 



(|201lD ); and (2) to p robe the strong gravitational field 
regime (jPsaltis 1 12008). However, although many models 
are available in the literature, the physica l origi n of kHz 
QPOs is still not known (e.g., iLin et al. I (|2011l) : iTorok I 
(|2009l) : lvan der Klisl (|2006D ). and hence we cannot yet use 



this promising feature as a tool. Many proposed mod- 
els for this t iming feature have a ttempted to explain the 
freque ncies dMiller et al. I Il998t IStella and Vietril Il998l 
1999; Lamb and Miller 20031: iKluzniak fc Abramowicz 
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Some of these models involve sev- 
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eral general relativistic frequencies at preferred radii and 
the neutron star spin frequency; in some cases beating 
and/or resonances among these frequencies. Other mod- 
els include the association of kHz QPOs with accretion 
through a non-axisy mmetric magnetic boundary lay er in 
the unstable regime omanova and KulkarnTI 12009): at- 
tempts to connect the kHz QPOs with the trapped two- 
armed nearly vertical oscillations in vertica l ly iso thermal 
disks with toroidal mag netic fields (|Kato 1 1201 ll ): attri- 
bution of the upper kHz QPO profile to the radial ex- 
tent of the kHz QPO emission region associated with 
the transitional layer a t the magnetosphere-disk bound- 
ary (jWang et al. 11201 ll) . etc. 

However, modeling only the frequencies gives an incom- 
plete picture; probing radiative transfer and further modu- 
lation processes is essential to constrain the existing mod- 
els for understanding the physical mechanism giving rise 
to kHz QPOs. For example, the models based on frequen- 
cies tentatively suggest some locations of kHz QPO origin. 
These locations could be at certain radii of the accretion 
disk, such as the innermo st-stable-circular-or bit radius, 
sonic point radius, etc. (see Ivan der Klisl (|2006| ) and refer- 
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ences therein). However, one needs to independently ver- 
ify these proposed locations. Since different energy spec- 
tral components originate from different locations, such 
as disk, boundary layer, corona, etc., a connection found 
between a kHz QPO property and a spectral component 
could provide this independent verification. Therefore, it 
is essential to study the energy dependence of kHz QPO 
properties. 

In this paper, we study the energy dependence of frac- 
tional rms amplitude of kHz QPOs. This property is 
a measure of QPO strength, and hence studying energy 
dependent kHz QPO property can be useful to probe 
which spectral component primarily contributes to this 
timing feature. This spectral connection of kHz QPO has 
been discussed by some authors. For example, iMendez I 
(2006) suggested that, while the kHz QPO frequencies 
are plausibly determined at the disk, this feature is 
modulated at the high-energy spectral component (e.g., 
corona, boundary layer, etc.). It has been reported by 
several authors that the kHz QPO fractional rms am- 
plitude increases with photon energy at lower ene rgies, 
and then plausibl y saturates (Ivan der Klisl (pOOl and 



references therein; iGilfanov et al. ( 20031 )1. Several au 



thors have theoretically computed the energy dependence 
of rms a mplitude of Com ptonizin g component variabil- 
ity (e .g .. ICabanac et alH (l2010t) : iGierlihski fc Zdziarskil 
(|2005l) ; IZdziarski et al.l (120051 ) : iLee and Miller I l|1998| )). 
For example, Cabanac et al. I (|2010l ) have demonstrated 
that oscillating hot thermal corona may give rise to an 
overall increase in rms variability with photon energy. 
However, the above mentioned saturation at higher en- 
ergies is not observationally established for many sources 
due to the lack of sufficient signal-to- noise ratios (S/Ns) 
at higher energies. For e xample, while th e high er energy 
data points of Fig. 10 of IGilfanov et al. I ((2003) are con- 
sistent with a flat fractional rms amplitude versus energy 
behavior, the errors of these data points are quite high. 

With this motivation, we try to measure the frac- 
tional rms amplitude versus energy behavior of lower kHz 
QPOs with improved S/N. We combine tens of lower kHz 
QPOs from the neutron st ar LMXB 4U 1728-34, in orde r 
to improve the S/N (e.g. , iFord and van der Kris I (11998 



Mc ndez and van der K lis (1999); Di Salvo et al. 



We choose lower kHz QPO, because this narrow QPO 
is more frequent and easier to detect than the relatively 
broad and weak upper kHz QPO. We find that the frac- 
tional rms amplitude systematically decreases with energy 
at higher energies. In order to understand this new find- 
ing, we compare the data with models involving a black- 
body+powerlaw energy spectrum. This comparison sug- 
gests that the observed behavior of fractional rms ampli- 
tude may be reproduced with a fluctuation of the black- 
body spectral component. In addition, we find that the 
quality factor (Q) of lower kHz QPO does not plausibly 
depend on energy. 

In § [2] and § O we describe our data analysis technique 
in detail and display our results, respectively. In § 21 we 
describe the models and discuss which ones are favored. 
Finally, in § [SJ we summarize our results and give impli- 
cations. 



2. DATA ANALYSIS 

In order to study the energy dependence of fractional 
rms amplitude associated with lower kHz QPOs, we 
retrieved all RXTE proportional counter array (PCA) 
pointed observational data of the source 4U 1728-34 for 
the period of April 13, 2000 to the end of December, 2009 
(« 272 ks of cleaned exposure time) corresponding to the 
PCA gain epoch 5. Note that gain, and hence the energy- 
channel conversion^, has changed during the entire RXTE 
lifetime. Therefore, we restrict our analysis to the data of a 
single epoch in order to make the energy dependence mea- 
surement more reliable. Although even within an epoch 
the energy-channel conversion evolves, such changes would 
be negligibly small to significantly affect the final conclu- 
sion of the analysis. We choose epoch 5, because it has the 
largest amount of data among all the epochs. We consider 
only the science event files with time resolution ~ 122/is 
and having continuous exposure time of at least 400 s for 
searching kHz QPOs. We do not apply any filtering based 
on energies, PCUs, Xenon layers, etc. However, following 
the usual practice, we filter the data based on time in or- 
der to remove thermonuclear X-ray bursts, data gaps, and 
observed intensity increase/decrease due to instrumental 
effects (especially due to start or stop of a PCU). 

It has been observed that all the kHz QPOs (which are 
strong and narrow) are confined to a small portion, i.e., 
low er banana, of the color- color diagram (CD) of 4U 1728- 
34 ((Di Salvo et al. 1120011 ). This suggests that they are 
of the same spectral origin. The strong int er-dependence 
found between Q-values and frequencies (|Barret et al. I 
2006!) of lower kHz QPOs from 4U 1728-34 suggests that 
all the QPOs considered by us are of same physical ori- 
gin. These justify the combination of many kHz QPOs (as 
mentioned in § [1} . 

We follow a few steps to probe the energy dependence of 
lower kHz QPO amplitudes. In the first step, we collect all 
the significant kHz QPOs from the event files (with at least 
400 s of data) using a blind sea rch. In order to do this, we 
compute a Leahy-normalized (|Leahv et al. I H9831) power 
density spectrum (PDS) using discrete Fourier transform 
DFTq from each 10 s of data for a given event file 
Ivan der Klis I [l989i)- The Nyquist frequency and fre- 
quency resolution of each such PDS (from 10 s of data ) 
are 2048 Hz and 0.1 Hz, respectively (|van der Klis II1989D . 
Then all (say, TV) the PDSs from a given event file are aver- 
aged, and we search for kHz QPOs in the frequency range 
of 400 — 1400 Hz in this mean PDS. In order to search 
effectively, several (say, W) adjacent frequency bins are 
combined. The noise powers in this PD S follow a \ 2 distr i- 
bution with 2 NW degrees of freedom (Ivan der Klis 1 11 989 h 
For a putative peak, we compute the single trial probabil- 
ity (q) of occurrence for a power greater than or equal to 
the peak power by chance from the noise power distribu- 
tion. Note that we consider the peak power, and do not fit 
the peak with a function (say, Lorentzian) for significance 
calculation. The probability q is multiplied with the num- 
ber of trials (Atrial) in order to obtain the final probability 
e (= gAt r i a i). We use A tr i a i = 20480, i.e., the original total 
number of frequency bins in a mean PDS. We consider a 
putative peak as a kHz QPO, if e < 4.65 x 10~ 4 (i.e., at 
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least 3.5cr significant). We detect tens of kHz QPOs from 
4U 1728-34 using the above procedure. 

In the second step, described in this paragraph, we iden- 
tify the lower kHz QPOs from the QPOs detected using 
the above mentioned criterion. If twin kHz QPO peaks 
appear in the mean PDS from an event file, the identifi- 
cation of the lower kHz QPO requires no further effort. If 
there is only one peak for an event file, then we consider it 
as a lower kHz QPO, if its Q > 40. Thi s is because, the Q 
vs frequency diagram of 4U 1728-34 (jBarret et al. II2006D 
clearly shows that only the lower kHz QPOs can have Q- 
values greater than 40, and all the upper kHz QPOs have 
smaller Q-values. In order to estimate the Q-values of 
kHz QPOs, we first u s e the shift-and - add te chnique (e.g., 
IMendez et al~l (|1998l ); iBarret et al~l 1)2006^ within each 
event file considering the entire PCA energy range. This 
technique mostly corrects for the centroid frequency drift. 
Note that such uncorrected drift makes the measured Q- 
value smaller than the actual value. After application of 
shift-and-add, we fit each PDS containing a detected kHz 
QPO with a 'constant + powerlaw + Lorentzian' model; 
constant takes care of the Poissonian white noise, pow- 
erlaw takes care of the low frequency red noise, and a 
Lorentzian describes the kHz QPO. We obtain the best-fit 
centroid frequency (v) and the best-fit full width at half 
maximum (FWHM) from the Lorentzian component, and 
compute the Q-value (i^/FWHM). Then we identify the 
lower kHz QPOs using the criterion Q > 40. From the 
entire « 272 ks of data, 40 lower kHz QPOs (each from 
one event file) are identified spanning 86.32 ks of exposure 
time, and are considered for further analysis. 

These 40 lower kHz QPOs are now to be combined in 
order to study the energy dependence of fractional rms 
amplitude with improved S/N. However, before doing this, 
we perform a few preliminary analyses in the third step, 
described in this paragraph. In this step, we still use the 
entire PCA energy range. We consider non-overlapping 
400-second segments from each of 40 event files with a 
lower kHz QPO. Now we collect only those 400-second 
segments, in each of which a lower kHz QPO peak with 
e < 2.7 x 10~ 3 (corresponding to 3ct) exists. To calculate 
the e for a given 400-second segment, we further divide the 
segment into 40 segments of 10-second intervals, calculate 
PDS for each 10-second interval, average these 40 PDSs, 
and use N = 40 and A/trial = 10000/ Note that this 
procedure, as well as the definitions of e, N, W and Atrial 
are given in the third paragraph of the current section. We 
find 173 number of 400-second segments with strong lower 
kHz QPO peaks using the above mentioned criterion on 
e. The mean of these peaks is v mcml = 807.8 Hz. We also 
calculate which is the separation between ^ mcan and 

the centroid frequency of an individual lower kHz QPO, 
for each 400-second segment. We use these 173 segments, 
fmean, and Vdifi values for further analysis. 

In the fourth step, described in this and the next three 
paragraphs, we compute lower kHz QPO fractional rms 
amplitudes for several chosen energy ranges after combin- 
ing the data of 173 number of 400-second segments (see the 
previous paragraph). We consider a set of PCA absolute 
channel ranges 5-8, 9-11, 12-13, 14-15, 16-17, 18-21, 20- 
25, 22-25, 22-29, 24-31, 26-33 and 26-49 corresponding to 



the energy ranges 2.06-3.68, 3.68-4.90, 4.90-5.71, 5.71- 
6.53, 6.53-7.35, 7.35-8.98, 8.17-10.63, 8.98-10.63, 8.98- 
12.28, 9.81-13.11, 10.63-13.93 and 10.63-20.62 keV, re- 
spectively, in order to study the energy dependence of 
fractional rms amplitude. We do not use proportional 
counter unit (PCU0) in this step, because the energy- 
channel conversion for PCU0 is somewhat different from 
that of the other four PCUs during epoch 5. For each en- 
ergy range, we compute a mean Leahy-normalized PDS (in 
the same way as described earlier in this section) for each 
of 173 number of 400-second segments. Then, we shift 
each of 173 lower kHz QPO peaks by v^is (the value for 
the entire PCA energy range; see the previous paragraph) 
to align them at i/ mcan , and add them together to obtain 
a grand average Leahy-normalized PDS for every energy 
range. 

Now the question is whether the centroid frequencies of 
the kHz QPOs change with energy. Such a change might 
affect our grand average Leahy-normalized PDS and fur- 
ther results. We cannot check it directly for individual kHz 
QPOs, because the statistics is not often adequate to de- 
tect an individual kHz QPO in a small energy range. So we 
check it in the following way. If we consider that centroid 
frequency of each of the individual kHz QPOs is energy 
dependent, then there are two extreme possibilities. (1) 
The first possibility: the centroid frequencies of all these 
kHz QPOs either increase or decrease at higher energy in 
a similar way. In this case, the i/ mcan value is also expected 
to either increase or decrease at higher energy, since we do 
not recalculate Vdifi separately for each energy range. But 
we find that the i^ mean value remains the same in all the 
energy bands (Figs. Q] and [2]). (2) The second possibility: 
the centroid frequencies of some kHz QPOs increase, and 
those of some other decrease with the increase of energy, 
in such a way that all the ^ moa n values at various energy 
ranges remain the same. In this case, although the cen- 
troid frequency would not change, the width of the merged 
kHz QPOs (i.e., after applying the shift-and-add method) 
is expected to systematically increase resulting in a lower 
and lower Q-value, at higher and higher energies. But we 
find that all the measured Q-values are quite close to each 
other (within the errors) and no systematic variation is 
observed. The fact that both v mean values and Q-values 
remain quite similar in all the energy ranges (Figs. [TJ [2] 
and suggests that the centroid frequencies of the kHz 
QPOs do not change with energy. 

The shift-and-add technique to obtain a grand average 
Leahy-normalized PDS (mentioned above) is a standard 
method to improve the S/N of kHz QPOs and to correc t 
their Q-value ( e. g.. IBarret et al. I (12006ft : IMendez I p)06); 
I Mendez et all (120011): IMendez and van der Klisl (I1999T) : 
IMendez et al. 
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(l2011tJ) b The grand average Leahy-normalized PDS for 
each energy range is then fitted with a 'constant + pow- 
erlaw + Lorentzian' model (as before), and the fractional 
rms amplitude and the Q-value are computed from the 
best-fit Lorentzian parameters (see the equation 4.10 of 
Ivan der Klis I (|1989f) for a general rms amplitude formula) . 

In the fifth step, we correct the above mentioned frac- 
tional rms amplitudes for the background levels. First, we 
compute the background count rate in each of the chosen 
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energy ranges using the appropriate channel range of the 
corresponding standard-2 data files. In order to do this, we 
use the PCA background model file for bright sources, and 
the 'FTOOLS' command 'pcabackest'. Then, we compute 
the background corrected fractional rms amplitudes (R c ) 
from the uncorrected fractional rms amplitudes (R uc ), the 
total count rates (I) and the background co unt rates (B) 
using the formula R c = R uc x 1/(1 — B) (M uno et al. I 
12002ft . 

Apart from the energy dependence of fractional rms am- 
plitude, we estimate how the Q-value of the lower kHz 
QPO depends on energy. This is computed from the above 
mentioned fitting of the grand average Leahy-normalized 
PDS for each energy range. The ratio of the best-fit cen- 
troid frequency to the best-fit FWHM of the Lorentzian 
model component gives the Q-value. Here we note that 
while the shift-and-add technique corrects for a large error 
in Q-value due to the centroid frequency drift, this tech- 
nique also introduces a small error in Q-value, because the 
centroid frequency is shifted without changing the FWHM. 

3. RESULTS 

The frequencies of the 173 detected lower kHz QPOs, 
which are used to study the energy-dependence of frac- 
tional rms amplitude (§0, span the range of 669.4 — 912.5 
Hz (mean (^ mca n) = 807.8 Hz; median = 810.2 Hz, and 
standard deviation = 57.8 Hz). In Figs. Q] and [2J the grand 
average Leahy-normalized PDSs (see § [5]) for 12 energy 
ranges are shown in 12 panels. Each panel is for the same 
exposure (69.2 ks), and in each panel the total count rate 
and the background count rate are mentioned. These fig- 
ures indicate that, while there are visibly prominent kHz 
QPO peaks for the lower energy bands, such a peak is not 
visible for the highest energy band. 

In Fig. [3l we plot background corrected fractional rms 
amplitude versus energy. The amplitude increases with en- 
ergy at lower energies, and systematically decreases (seen 
from the overlapping energy ranges) above « 10 keV. Note 
that, if for an energy range the lower kHz QPO peak 
is not clearly visible (see Fig. [3J panels with blue model 
curves), we compute an upper limit to the fractional rms 
amplitude at i-Wan (see § [2]) instead of estimating the am- 
plitude from the best-fit Lorentzian. Our estimated ler 
and 3(7 upper limits of fractional rms amplitude in the en- 
ergy range 10.63-20.62 keV are 2.37% and 3.66%, respec- 
tively. If we compare the fractional rms amplitude value 
(= 14.38% ± 1.71%) in the 7.35-8.98 keV range with the 
3(7 upper limit (= 3.66%) in the 10.63-20.62 keV range, we 
find a 6.3cr significant drop in the amplitude value. This 
quantifies the decrease of the lower kHz QPO fractional 
rms amplitude at higher energies (see Fig. [3]). 

The initial increasing part of the fractional rms ampli- 
tude is fitted with a linear and a powerlaw model, the 
best-fit curves of which are very similar to each other 
(see Fig. 0]). This figure, with extrapolated model curves, 
shows that the increasing trend of lower kHz QPO ampli- 
tude at lower energies cannot plausibly explain the upper 
limit in the 10.63-20.62 keV range. For example, the ex- 
trapolated fractional rms amplitude value at 10.63-20.62 
keV from the linear model is 30.08 ± 6.80, which implies 
a 3.89cr significant higher value compared to the 3cr upper 
limit of the observed fractional rms amplitude. 

While the fractional rms amplitude changes with energy, 



Q-value does not (Fig. [31 lower panel). As mentioned in 
§ [2j while the shift-and-add technique corrects for a usually 
large error in Q-value due to the evolution of the centroid 
frequency, this technique introduces a small error. This 
is because, while a kHz QPO peak is shifted to a new 
frequency (u meaa = 807.8 Hz), the FWHM of this peak 
remains the same. Given that the standard deviation of 
the centroid frequency distribution of our detected lower 
kHz QPOs is 57.8 Hz and tWan = 807.8 Hz, this error can 
be estimated to be ~ 7% of the typical value of Q ~ 80. 
Therefore, this error is small compared to the la error 
(which is typically ~ 25% — 45%) from Lorentzian fitting 
of the QPO peaks. Therefore, the error due to frequency 
shift does not change our conclusion regarding the energy 
independence of Q-values, and hence we do not make an 
attempt to correct it. 

4. COMPARISON WITH MODELS 

This paper shows a prominent and systematic decrease 
of the lower kHz QPO fractional rms amplitude at higher 
energies, for the first time to the best of our knowledge. 
We, therefore, try to find out if such a decrease is at all the- 
oretical expected. For example, could certain parameter 
values of the usual energy spectra of neutron star LMXBs 
explain it qualitatively? Before testing this, let us briefly 
mention the X-ray components of these sources, as we un- 
derstand currently. Neutron star LMXBs are believed to 
have two primary X-ray emitting regions, an accretion disk 
and a boundary layer. Both these regions are expected to 
be optically thick, and therefore to emit blackbody radi- 
ation. Furthermore, one (or more) of these components 
may be covered with a corona (coronae). Such a corona 
may reprocess (Comptonize) the blackbody photons. The 
amount of reprocessing depends on the optical depth of the 
corona, as well as the extent of coverage (e.g., full versus 
partial). As a result, the observed spectrum either from 
the disk and/or the boundary layer can be a blackbody, or 
Comptonized, or a sum of both. Observationally, however, 
no spectral model usually uniquely describes the data. 
Moreover, it is not usually clear, where various compo- 
nents of a given model originate from. The refore, although 
many wo rks over a few deca des (e.g. , IMitsuda et al. ' 



IMitsuda and Dotani I (U9 89) ; 



Christian and Swank 



( 2003) ; Maitra and Bailvnl (120041); iGilfanov et al. 



White et al . 

Church and Baluririska-Church I (120011)" 
1997); Mac carone and Coppi 



2003 ): lOlive et al l (l200ak IWiinands I (120011): iBarret 



j200ll);lLin et al. Tfl2007ll201C|) : lMukheriee and Bhattacharyya I 
(|2011af) ) support the general understanding mentioned 
above, the details are still unknown. 

In order to have a basic understanding (i.e., as much 
as possible independent of detailed models) of the lower 
kHz QPO fractional rms amplitude versus energy behav- 
ior found by us, we choose the simplest two-component 
spectral model, which is based on the above mentioned 
description of X-ray components. This is a black- 
body+powerlaw model, in which po werlaw usually repre - 
sent s the Comptonized componen t. iTarana et ail ([201 ll ) 
and lSeifina and Titarchukl (|201ll ) have successfully fitted 
the continuum spectra from 4U 1728-34 with a black- 
body+Comptonization model. In order to check whether 
a blackbody+powerlaw works, we fit the continuum spec- 
tra of a few data files (having prominent kHz QPOs) with 
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this model, as well as with a blackbody-f-Comptonization 
(bbodyrad+compTT in XSPEC) model. We find that, not 
only both the models give acceptable fits, but the power- 
law component represents the Comptonization component 
reasonably well, as verified from the respective fluxes and 
spectral component curves. This indicates that a black- 
body+powerlaw model is appropriate for the 4U 1728-34 
spectra. 

Before describing our model, here we mention and dis- 
cuss the assumptions involved in our modeling and inter- 
pretation, some of them recapitulated from the data analy- 
sis side and some other extra assumptions from the model 
side. (1) We assume that the spread in energy spectra 
(e.g., soft colors and hard colors in a color-color diagram 
(CD): Ivan der Klisl ()2006f )) and in intensity from our data 
files containing lower kHz QPOs do not significantly affect 
the robustness of the energy behavior of fractional rms 
amplitude reported by us (Fig. [3]). This is justified be- 
cause of the following reasons, (a) All the QPOs (which 
are strong and narrow) considered by us are confined to a 
small po rtion, i.e., lower banana, of the CD of 4U 1728- 
34 (|Di Salvo et al. II2001D . as mentioned in §H (b) Given 
the small fractional rms amplitudes (~ 5 — 14%) of lower 
kHz QPOs in the 2-16 keV energy range, which is used to 
compute the CD and HID, the entire X-ray energy spec- 
trum does not contribute to these QPOs. Hence a small 
spread in CD does not necessarily indicate different phys- 
ical origin of different lower kHz QPOs. (c) The strong 
inter-dependence found between Q-values and frequencies 
of lower kHz QPO s from 4U 1728-34 ( along with several 
other atoll sources: iBarret et al. I ()2006[ )) suggests that all 
the QPOs considered by us are of same physical origin. 
These support the robustness of our reported rms-energy 
behavior. (2) We assume that the energy dependence of 
lower kHz QPOs does not change with frequency. This 
may be a reasonable assumption, because these QPOs with 
different frequencies are likely to have originated from the 
same physical process. (3) Our model considers the si- 
nusoidal fluctuation of a single parameter of the above 
mentioned blackbody+powcrlaw model. The sinusoidal 
fluctuation is somewhat justified, because according to 
some models, the original oscillations are sinusoidal, and 
the broadening happe ns because o f a dec oherence mech- 
anism (say, damping; Ivan der Klisl (|2006l) and references 
therein). H ere we note tha t a sin usoidal signal was also 
assumed by I Lee and Millerl (fl99 8). We further note that, 
even if the intrinsic oscillations of lower kHz QPOs are 
a sum of many sinusoids originated separately, our mod- 
eling will still be useful to connect a spectral component 
to these QPOs. (4) Our model does not involve any re- 
sponse matrix. Although application of a response matrix 
may slightly change the model rms-energy curves quanti- 
tatively, the qualitative nature of these curves should not 
change. 

Now we describe our model. As previously mentioned, 
the total time averaged flux (say, < S(E,t) >; E: pho- 
ton energy, t: time) of the source can be fitted well with 
a combination of blackbody (Sbb) and powerlaw (Spl). 
Each of blackbody and powerlaw has two parameters: nor- 
malization (Abb) and temperature (Tbb) for the former, 
and normalization (Apl) and photon index («pl) for the 
latter. We try to investigate if sinusoidal fluctuation in 



one of these parameters (while keeping the other three 
parameters non-fluctuating) can reproduce the observed 
fractional rms amplitude versus energy behavior. Here as 
an example, we consider the case of fluctuating blackbody 
temperature. 

N H x S(E,t) = A H x [S BB (E,N BB .T BB (t)) 

+ S PL (E,Np L ,a PL )} 
where, T B b(£) = T B bo + T BB1 Sin(ujt) 
A H x < S(E, t) > = A H x [< S BB (E, N BB , (T BB0 

+ T Bm Sin(ut))) > +S PL (E,Np L , apL )} 
A H x a s (= <7 NS ) = A H x V(< S(E,t) 2 >-< S(E,t) > 2 ) 

fracRMS = r — 

A H x < S{E,t) > 

V(< S(E,t) 2 >-<S(E,t) > 2 ) 
< S(E,t) > 

which is independent of Ah (neutral hydrogen column den- 
sity). Here, fracRMS is our model fractional rms ampli- 
tude. Similarly, we compute the model fractional rms am- 
plitude in each case when any one of the other three pa- 
rameters varies sinusoidally in time. 

We note that the blackbody flux and the powerlaw 
flux arc nonlinear functions of Tbb and cupl, respectively. 
Therefore, if Tbb (or apL) fluctuates sinusoidally (as we 
assume) , the variation of the blackbody flux (or powerlaw 
flux) will not be strictly sinusoidal, and will contain higher 
harmonics. Here we show that this flux variation is sinu- 
soidal (which we assume in our modeling) when Tbb (or 
«pl) fluctuation is small (e.g., Tbbi *C Tbbo, for the vary- 
ing blackbody temperature), and is therefore considered 
up to the first order: 

< S{E,t) > = < S BB (E,N BB , (T BB0 + T BB1 Sin(ut))) > 
+ < S Ph {E,Np h ,a PL ) > 
— < S BB (E,N BB ,T BB0 ) > 
+ < S B m{E, N BB ,T BB0 ,T BB i)Sin(u}t) > 

+ 0[^I] 2 + S PL (£,Ap L ,ap L ) 
Jbbo 

= S BB (E, A B b,T B bo) + Spl(E, A PL ,ap L ) 

r T BB i 



(which is correct up to first order in [ 



bbo 



The last step is justified because the observations are sev- 
eral orders of magnitudes longer than the oscillation pe- 
riod T = 2-k/uj. This last step shows that, < S(E,t) > 
= S(E); which is simply the sum of XSPEC model compo- 
nents (bbodyrad+powerlaw) used to fit the time averaged 
spectrum. Finally, note that this first order approxima- 
tion is not required when normalization of blackbody or 
normalization of powerlaw varies. 

Considering the above procedure, our model fractional 
rms amplitude has five parameters: Abb, Tbb, Apl, Q:pl 
and A, the last one is the fractional peak amplitude of 
one of the first four parameters, which is fluctuating (e.g., 
y BB1 in the case discussed above for the varying blackbody 
temperature). While computing the models, we consider 
following ranges of our parameters: Abb =0.5— 50, Tbb = 
0.3-3.0 keV, A PL = 0.5-50 photons cm" 2 s" 1 keV" 1 (at 
1 keV), a PL = 1.5-3.5 and A = 0.001-0.9. We note that 



6 



Mukherjee and Bhattacharyya 



the formula for blackbody photon count rate per cm 2 per 
keV is (A^BBxl.0344xlO- 3 x£; 2 )/(exp(£;/TBB)-l), while 
that for powerlaw photon count rate per cm 2 per keV is 
-/VpL-E~ apL , with energy E in keV. The best-fit parameter 
values from our fitting of a few 4U 1728-34 spectra (see 
earlier in this section), and fr om lTarana et al. I (| 2 1 lh and 
ISeifina and Titarchuk I (|201lh are within the above ranges. 

For a set of five model-parameter values, we numerically 
compute the fractional rms amplitude, as outlined above, 
systematically in each of the observed energy ranges for 
the sinusoidal oscillation in time of each of A^bb, Tbb, ^Vpl 
and cupL- For each of these four cases, we produce many 
'fractional rms amplitude versus energy' model curves for 
various sets of parameter values, in order to thoroughly 
explore if some curves for a given case can reasonably de- 
scribe the data. We show some examples of these model 
curves in Figs.[51[IO[Z]and[5] In each of these figures, for the 
variation of a single model parameter, we show the exam- 
ple curves for various sets of parameter values (mentioned 
in the figure-caption) . Therefore, these figures give an idea 
about the nature of energy-dependence of fractional rms 
amplitudes for various parameter values for the sinusoidal 
oscillation in time of a given parameter. The blue dashed 
curve in each figure appears to be the closest to the data 
among our computed model curves. These figures show 
that the sinusoidal oscillation of a blackbody parameter 
describes the data well, and hence the observed energy- 
dependence of fractional rms amplitude is expected. We 
note that testing this was the main aim of modeling re- 
ported in this paper, as mentioned in the first paragraph 
of this section. Figs. El El [7] and [8] also suggest that the os- 
cillation of a single blackbody parameter (temperature or 
normalization) describes the data better than the variation 
of a single powerlaw parameter (index or normalization). 
We further note that the small value of A corresponding 
to the blue dashed curve of Fig. [5] satisfies the assump- 
tion of small fluctuation of the blackbody temperature, as 
mentioned earlier in this section. 

5. DISCUSSION 

In this paper, we report the decrease of fractional rms 
amplitude of lower kHz QPOs at higher energies, for the 
first time to the best of our knowledge. As the count rate 
of neutron star LMXBs decreases at higher energies, it is 
difficult to reliably measure the fractional rms amplitude 
with RXTE PCA (the only instrument capable of detect- 
ing kHz QPOs) due to the lack of sufficient S/N. Some- 
times because of this, the previous studies, to the best 
of our knowledge, have concluded that the fractional rms 
amplitude either goes on increasing with energies, or first 
increases and then saturates. In order to improve the S/N, 
we combine tens of lower kHz QPOs from 4U 1728-34, and 
also use the standard shift-and-add technique (§H]). The 
net count rate of our highest energy range (10.63 — 20.62 
keV) is also larger than that in some other energy ranges. 
Since the 10.63 — 20.62 keV range does not have any clear 
kHz QPO peak, the decrease of fractional rms amplitude 



at higher energies should be real, and not a result of low 
S/N (see Figs. Q] and [2]). Note that the small uncertainties 
in the channel-energy conversion cannot make the lower 
kHz QPO disappear at higher energies, and hence cannot 
produce the decrease. Moreover, the gradual and system- 
atic change of fractional rms amplitude at higher energies 
for overlapping energy ranges (Figs. [HOE]) strongly sup- 
ports the decrease of this amplitude of lower kHz QPOs. 

The energy dependence of fractional rms amplitude can 
be very useful to probe the physical origin of kHz QPOs 
(see §H|), and hence our finding is important. Therefore, 
we use a two-component (blackbody+powerlaw) spectral 
model to gain insights (see §2] for details). We explore the 
shape and fractional rms amplitude values of many model 
curves, examples of which have been shown in Figs. [5j O 
and[H] These curves suggest that (1) the reported energy 
behavior of lower kHz QPO fractional rms amplitude can 
be explained with simple models; and (2) it is likely that 
the fluctuation in a blackbody-like component primarily 
causes the lower kHz QPOs. However, we consider the 
fluctuation of only one model parameter at a time (§ 3]), 
and hence we cannot study the effects of a simultaneous 
fluctuation of the two powerlaw parameters. We also note 
that the non-dependence of Q-value on energy (§ ^ may 
imply that the lifetime of the oscillations is independent 
of photon energies. This may suggest that the lower kHz 
QPOs in various energy ranges originate from the same lo- 
cation, and hence from the same spectral component. But 
we cannot constrain the location of either the blackbody or 
the powerlaw component from our study. The blackbody 
could originate eithe r from the boundary layer or the disk 
(|Tarana et al. 11201 ll) . Similarly, the corona, which plausi- 
bly gives rise to the powerlaw, could be centrally located 
near the neutron star or a cover on the accretion disk. 

With our new finding, an important question to ask 
would be whether the fractional rms amplitude of lower 
kHz QPOs decreases with energy (at higher energies) for 
other neutron star LMXBs as well. Gilf anov et al. 1 |2003) 
reported that 4U 1608-52 did not show any decrease in 
fractional rms amplitudes at energies as high as 20 keV, 
although at the high energy end the error in fractional rms 
amplitude was quite high. A larger (than PCA) area at 
higher energies may be required (1) to detect a plausible 
decrease which happens at higher energies than that for 
4U 1728-34; (2) and to reduce the errors in fractional rms 
amplitudes, while increasing the number of data points, in 
order to perform a better modeling. The upcoming space 
missions (e.g., ASTROSAT) having sufficient time resolu- 
tion and much higher (than PCA) area in 20 — 50 keV, will 
be useful for this purpose. 

We thank Manoneeta Chakraborty and H. M. Antia for 
technical discussions. A.M. acknowledges Ranjeev Misra 
for a discussion on the models and thanks M. Coleman 
Miller for useful general discussions on kHz QPOs. We also 
thank an anonymous referee for very constructive com- 
ments, which significantly improved the paper. 
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Fig. 1. — Leahy-normalized power density spectra (PDS; with lcr errors) of 4U 1728-34 after using the shift-and-add technique for each of 
the energy ranges 2.06-3.68 keV (panel a), 3.68-4.90 keV (panel 6), 4.90-5.71 keV (panel c), 5.71-6.53 keV (panel d), 6.53-7.35 keV (panel 
e) and 7.35-8.98 keV (panel /; see §E}. In each panel, the average total count rate (i.e. without background subtraction) used to compute 
the PDS is denoted with 'total', the average background count rate is denoted with 'bkg', and the PCA absolute channel range is denoted 
with 'Ch'. Furthermore, the frequency resolution and exposure in each panel are 3.2 Hz and 69.2 ks, respectively. A lower kHz QPO peak is 
clear lv seen in each Danel. The best-fit 'constant + Dowcrlaw + Lorcntzian' model for each enerev ranee is shown with a red curve (§\2\i. The 
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Fig. 2.— Same as Fig.[TJ but for the energy ranges 8.17-10.63 keV (panel a), 8.98-10.63 keV (panel 6), 8.98-12.28 keV (panel c), 9.81-13.11 
keV (panel d), 10.63—13.93 keV (panel e) and 10.63-20.62 keV(panel /; see §[2j. In the panels d, e and /, the lower kHz QPO peak is not 
clearly seen, and hence the best-fit model curves in these panels are marked with a different color (blue). For each of these three energy 
ranges, an upper limit to the lower kHz QPO fractional rms amplitude is computed (see §[3j. 
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Fig. 3. — Energy dependence of the background corrected fractional rms amplitude and the Q-value of the lower kHz QPO from 4U 1728—34. 
In each panel, solid horizontal lines denote a set of energy ranges which are mutually non-overlapping, and energy ranges denoted with dotted 
horizontal lines overlap with some ranges denoted with solid lines. Upper panel: rms amplitude (with la error or la upper limit) versus 
energy with ranges 2.06-3.68 keV (Fig. Q] panel a), 3.68-4.90 keV (Fig. [TJ panel b), 4.90-5.71 keV (Fig. Q] panel c), 5.71-6.53 keV (Fig. Q] 
panel d), 6.53-7.35 keV (Fig. Q] panel e), 7.35-8.98 keV (Fig. Q] panel f), 8.17-10.63 keV (Fig. EJ panel a), 8.98-10.63 keV (Fig. El panel 
b), 8.98-12.28 keV (Fig. El panel c), 9.81-13.11 keV (Fig. El panel d) and 10.63-20.62 keV (Fig. El panel f). This panel shows that the rms 
amplitude clearly and gradually decreases at higher energies. Lower panel: Q-value (corrected for frequency drift with shift-and-add technique 
for 400 s segments; with la error) versus energy with ranges 2.06-3.68 keV (Fig. [l] panel a), 3.68-4.90 keV (Fig.[l] panel b), 4.90-5.71 keV 
(Fig. [I] panel c), 5.71-6.53 keV (Fig. Q] panel d), 6.53-7.35 keV (Fig. □ panel e), 7.35-8.98 keV (Fig. Q] panel f), 8.17-10.63 keV (Fig. El 
panel a), 8.98-10.63 keV (Fig. El panel b) and 8.98-12.28 keV (Fig. El panel c). This panel shows that the Q-value does not significantly 
depend on energy (see §E|and[3]|. 
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Fig. 4. — Energy dependence of the background corrected fractional rms amplitude of the lower kHz QPO from 4U 1728-34. The first six 
data points (with la errors), which are same as the first six data points of the upper panel of Fig. [3] are fitted with a linear (solid blue line) 
and a powerlaw model (dotted red line). Both the best-fit models (extrapolated) are shown. The last data point (for the same energy range 
as that of the last data point of the upper panel of Fig. [3j gives the 3cr upper limit of fractional rms amplitude. This figure shows that the 
increasing trend of amplitude at lower energies cannot explain the last data point. 
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Fig. 5. — The lower kHz QPO fractional (in %) rms amplitude vs. energy. The red points are the same non-overlapping data points of the 
upper panel of Fig. [3] (§ [2] and [3)1 . The curves are models for sinusoidal fluctuation of the blackbody temperature (see §HJ. The parameter 
(^Vpl, c*PL> ^BBi Tbb and A) values for various model curves are as follows in the same sequence (for units see § [4j . 1: [0.5, 3.5, 50.0, 0.89, 
0.034]; 2: [0.5, 3.5, 0.5, 0.88, 0.034]; 3: [0.5, 1.5, 50.0, 0.88, 0.034]; 4: [0.5, 3.5, 50.0, 3.0, 0.034]; 5: [0.5, 3.5, 50.0, 1.5, 0.034]; 6: [0.5, 3.5, 
50.0, 0.5, 0.034]; 7: [0.5, 3.5, 50.0, 0.88, 0.1]; 8: [0.5, 3.5, 50.0, 0.88, 0.01]. The blue dashed curve (curve 1) appears to be the closest to 
the data among our computed model curves for blackbody temperature fluctuation (see §|3). This figure shows that blackbody temperature 
fluctuation could reproduce the observed energy dependence of the lower kHz QPO fractional rms amplitude (§[4]|. 
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Fig. 6. — Similar to Fig. [5] but the curves are models for sinusoidal fluctuation of the blackbody normalization (see §|4jl. The parameter 
(■^PLi cipL, ^VbBi Tbb an( i -A) values for various model curves are as follows in the same sequence (for units see §@J. 1: [3.69, 2.79, 10.37, 
1.42, 0.9]; 2: [0.5, 3.5, 13.8, 1.11, 0.2]; 3: [0.5, 3.5, 13.8, 1.11, 0.5]; 4: [0.5, 3.5, 12.4, 3.0, 0.1]; 5: [0.5, 3.5, 50.0, 3.0, 0.05]; 6: [10.0, 1.5, 10.0, 
2.0, 0.1]; 7: [10.0, 3.5, 10.0, 1.0, 0.1]; 8: [0.5, 2.0, 10.0, 1.5, 0.1]; 9: [50.0, 3.5, 50.0, 3.0, 0.25]. The blue dashed curve (curve 1) appears to be 
the closest to the data among our computed model curves for blackbody normalization fluctuation (see §[2). This figure shows that blackbody 
normalization fluctuation could reproduce the observed energy dependence of the lower kHz QPO fractional rms amplitude, although for a 
high value of A (§0). 
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Fig. 7. — Similar to Fig. [5] but the curves are models for sinusoidal fluctuation of the powcrlaw index (sec §2}. The parameter (-/Vpl, apL, 
^BBi ^bb an d A) values for various model curves are as follows in the same sequence (for units see §0}. 1: [18.55, 3.5, 2.07, 3.0, 0.027]; 2: 
[45.0, 3.5, 5.0, 2.5, 0.03]; 3: [45.0, 3.5, 5.0, 2.0, 0.03]; 4: [45.0, 3.5, 5.0, 1.0, 0.03]; 5: [45.0, 3.5, 50.0, 3.0, 0.03]; 6: [45.0, 3.5, 50.0, 2.0, 0.03]; 7: 
[10.0, 3.5, 5.0, 3.0, 0.03]; 8: [10.0, 2.0, 10.0, 1.5, 0.001]; 9: [10.0, 2.0, 10.0, 1.5, 0.1]. The blue dashed curve (curve 1) appears to be the closest 
to the data among our computed model curves for powerlaw index fluctuation (see §|4j. This figure shows that powerlaw index fluctuation 
plausibly could not reproduce the observed energy dependence of the lower kHz QPO fractional rms amplitude (see §|4{. 
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Fig. 8. — Similar to Fig. [5] but the curves are models for sinusoidal fluctuation of the powcrlaw normalization (see §|4)l. The parameter 
(TVpL, «pli ^BBi Tbb an d ^4) values for various model curves are as follows in the same sequence (for units see §U)l. 1: [0.5, 3.5, 50.0, 0.52, 
0.117]; 2: [10.0, 2.0, 10.0, 1.5, 0.02]; 3: [10.0, 3.5, 10.0, 1.5, 0.2]; 4: [10.0, 2.0, 50.0, 1.5, 0.2]; 5: [50.0, 3.5, 50.0, 3.0, 0.2]; 6: [0.5, 3.5, 0.5, 
2.0, 0.2]; 7: [0.5, 3.5, 50.0, 1.5, 0.12]; 8: [0.5, 3.5, 50.0, 0.5, 0.3]; 9: [0.5, 3.5, 50.0, 1.5, 0.3]. The blue dashed curve (curve 1) appears to be 
the closest to the data among our computed model curves for powerlaw normalization fluctuation (see ? [4l . This figure shows that powcrlaw 
normalization fluctuation plausibly could not reproduce the observed energy dependence of the lower kHz QPO fractional rms amplitude (see 

§U>- 



